A = {"/tmp/tmp.cfg", "00011"} ;
Yfut = [2071 1 1; 2100 12 31] ;
sfx = "z" ; sfun = @mean ; rot = 90 ;

global PAR PRJ isANA ADIR CDIR SDIR isCUR

addpath([pwd "/fun"]) ;
addpath([pwd "/util"]) ;
Lptr = Lpdd = Lcal = Lout = Lplt = true ;

for i = 1:numel(A)
   if exist(A{i}, "file")
      A{i} = ["CFG = \"" A{i} "\";"] ;
   elseif regexp(A{i}, "^[01]{5}$")
      [Lptr Lpdd Lcal Lout Lplt] = deal(cellfun(@(x) x == "1", num2cell(A{i}), "UniformOutput", false){:}) ;
      A{i} = {} ;
   endif
endfor

init(A{:}) ;                     # initialization

pkg load general miscellaneous

ldata = fullfile("data", PRJ, PAR.pdd) ;

printf("ptr_ana: <-- %s\n", ofile = fullfile(ldata, "pdd.ob")) ;
load(ofile) ;
printf("ptr_ana: <-- %s\n", cfile = fullfile(ldata, "clim.ob")) ;
load(cfile) ;
printf("ptr_ana: <-- %s\n", PAR.mfile) ;
load(PAR.mfile) ;
if strcmp(PAR.mod, "xds")
   S = E ;
else
   S = R ;
endif

vars = pdd.vars ; VAR = sunique(vars) ;

### simulation
printf("ptr_ana: <-- %s\n", afile = fullfile(ADIR, [strrep(PAR.ana, "/", "_") ".pc.ob"])) ;
[xds X] = dsc(afile, clim, S) ;
[ana.id ana.x] = annstat(xds.id, xds.(sfx), sfun) ;

printf("ptr_ana: <-- %s\n", vfile = fullfile(SDIR, [strrep(PAR.ana, "/", "_") ".pc.ob"])) ;
[xds X] = dsc(vfile, clim, S) ;
fut.PC = selper(X, sdate(X.id, Yfut(1,:), Yfut(2,:))) ;
[fut.id fut.x] = annstat(xds.id, xds.(sfx), sfun) ;

printf("ptr_ana: <-- %s\n", v0file = fullfile(CDIR, [strrep(PAR.ana, "/", "_") ".pc.ob"])) ;
[xds X] = dsc(v0file, clim, S) ;
cur.PC = selper(X, sdate(X.id, PAR.cper(1,:), PAR.cper(2,:))) ;
[cur.id cur.x] = annstat(xds.id, xds.(sfx), sfun) ;

Xs = mean(fut.PC.x) - mean(cur.PC.x) ;

## [In Ix] = ext(xds.(sfx)) ; 
## Xs = X.x(Ix(1),:) ;

if PAR.pcr
   load(strrep(PAR.mfile, "mod", "pcr")) ;
   Xs = Xs * eof.E ; S = mult(eof.E', S) ;
   J = find(abs(S(:,PAR.validx,1)) < eps) ;
   tl = [repmat("EOF ", length(J), 1) num2str(J)] ;
else
   Xj = Xjfun(cur.PC.J) ; t = find(Xj == 1) ; tl = cur.PC.ptr(cur.PC.J)(t) ;
endif

set(0, "defaulttextfontsize", 8, "defaultaxesfontsize", 8) ;
figure(1) ;

PAR.validx = 4 ;

L = isempty(PAR.validx) ;
n = L*pdd.nc + !L*numel(PAR.validx) ;
ylim = [-7 7] ; nc = min(3, n) ; nr = ceil(n/nc) ; k = 0 ;
for j=1:pdd.nc

   if !L && j != PAR.validx, continue ; endif
   k++ ;
   subplot(nr, nc, k) ;
   plot(cumsum(Xs .* S(:,j)'), "color", "black") ;
   plot(nancumsum(Xs .* S(:,j)'), "color", "blue") ;
   set(gca, "XTick", t, "XTickLabel", {tl{:}}, "xgrid", "on", "ygrid", "on") ;
   if L && strcmp(vars{j}, "P")
      set(gca, "ylim", ylim) ;
   endif
   legend(sprintf("%s, %s", pdd.ids{j}, vars{j})) ;
   legend("location", "northwest") ;
   rotateticklabel(gca, rot) ;

endfor

figure(2) ;
clf ; j = PAR.validx ;
plot(ana.id(:,1), ana.x(:,j), "b", cur.id(:,1), cur.x(:,j), "g", fut.id(:,1), fut.x(:,j), "r") ;

figure(3) ;
clf ;
subplot(311) ; plot(S(:,j)) ; title("S") ;
set(gca, "XTick", t, "xgrid", "on", "ygrid", "on") ;
subplot(312) ; plot(Xs') ; title("Xs") ;
set(gca, "XTick", t, "xgrid", "on", "ygrid", "on") ;
subplot(313) ; plot(Xs' .* S(:,j)) ; title("Xs*S") ;
set(gca, "XTick", t, "XTickLabel", {tl{:}}, "xgrid", "on", "ygrid", "on") ;
rotateticklabel(gca, rot) ;
exit

clf ;
iptr = 1 ; j = Xj(t(iptr)) ;
load(fullfile(gdata, PAR.ana, ANASFX, "eof", [GVAR{1,iptr} ".ob"])) ;
load(fullfile(gdata, PAR.sim{2}, "ptrz.ob")) ;
pc.x = ptr.(GVAR{2,iptr}).z * eof.E ;
[pc2.id pc2.x] = annstat(ptr.(GVAR{2,iptr}).id, pc.x, sfun) ;
load(fullfile(gdata, PAR.sim{1}, "ptrz.ob")) ;
pc.x = ptr.(GVAR{2,iptr}).z * eof.E ;
[pc.id pc.x] = annstat(ptr.(GVAR{2,iptr}).id, pc.x, sfun) ;
plot(pc2.id(:,1), pc2.x(:,j), "g", pc.id(:,1), pc.x(:,j), "r") ;

## plot pattern
addpath ~/matlab/m_map ~/matlab/m_map/private

figure ;
clf
m_proj('albers','long',PAR.lon,'lat',PAR.lat, "par", affin(PAR.lat, 0.5));
m_coast('patch',[1 .85 .7]);
[cs, h] = m_elev('contour',[100:200:1000]) ;
m_grid('box','fancy','tickdir','in');
colormap(flipud(copper));

[lon,lat]=meshgrid(eof.lon, eof.lat);
eof.E = reshape(eof.E, length(eof.lat), length(eof.lon), []) ;
hold on
[cs,h]=m_contour(lon, lat, eof.E(:,:,1));
clabel(cs, h, 'fontsize', 8);
xlabel('Simulated something else');


%% Nice looking data
[lon,lat]=meshgrid([-136:2:-114],[36:2:54]);
u=sin(lat/6);
v=sin(lon/6);

m_proj('oblique','lat',[56 30],'lon',[-132 -120],'aspect',.8);

subplot(121);
m_coast('patch',[.9 .9 .9],'edgecolor','none');
m_grid('tickdir','out','yaxislocation','right',...
       'xaxislocation','top','xlabeldir','end','ticklen',.02);
hold on;
m_quiver(lon,lat,u,v);
xlabel('Simulated surface winds');

subplot(122);
m_coast('patch',[.9 .9 .9],'edgecolor','none');
m_grid('tickdir','out','yticklabels',[],...
       'xticklabels',[],'linestyle','none','ticklen',.02);
hold on;
[cs,h]=m_contour(lon,lat,sqrt(u.*u+v.*v));
clabel(cs,h,'fontsize',8);
xlabel('Simulated something else');


clf
m_proj('lambert','long',[-120 -100],'lat',[30 60]);
m_coast('patch',[1 .85 .7]);
m_elev('contour',[500:500:6000]);
m_grid('box','fancy','tickdir','in');
colormap(flipud(copper));
